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ABSTRACT 



Subsonic and transonic steady and unsteady flowfields over airfoils axe inves- 
tigated with the numerical solution of the governing equations. This study aims to 
enhance the performance of rotary wing and fixed wing aircraft by better understand- 
ing and by talcing advantage of unsteady phenomena such as dynamic lift. In the past 
few years many advances have been made in algorithm development for the numerical 
solution of the Euler and the Navier Stokes equations. In this study, these new zonal 
techniques axe applied. A zonal approach is more computationally efficient in solving 
the governing equations than previous approaches, and has certain advantages over 
the standard single moving grid approach. The zonal grid consists of two grids, one 
being the inner grid which is fixed to the airfoil, and the other being the outer grid 
which extends to the fax field or to a specified outer boundary. The inner grid is 
allowed to rotate with the body, while the outer grid remains fixed. The thin-layer 
Navier-Stokes equations axe solved for the inner grid, while the Euler equations are 
solved for the outer grid. Communication between the two grids is accomplished by 
interpolating the flow quantities at the zonal interface. Solutions axe obtained for 
flows at fixed angles of incidence, and for unsteady flows over pitching and oscillating 
airfoils. The computed results are in good agreement with available experimental 
data. 
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I. INTRODUCTION 



A. BACKGROUND 

Investigation of steady and unsteady flowfields over airfoils is an active area 
of numerical and experimental research. Unsteady pitch-up motion of airfoils alters 
significantly the aerodynamic characteristics of lifting surfaces. Pitch-up motion of 
an airfoil produces lift augmentation and delay of stall to higher angles of incidence 
compared to airfoils held at a steady incidence. Understanding the mechanisms that 
cause the dynamic lift development is a subject of interest in both theoretical and 
applied research. 

A comprehensive explanation of the dynamic lift phenomenon is given by Tyler 
and Leishman [Ref. 1]. To briefly summarize, unsteady airfoil motion tends to 
maintain high lift to higher angles-of-attack, because the unsteady flow causes a 
time delay in the build-up of the lift force and the adverse pressure gradient. The 
unsteady motion also gives the airfoil a virtual camber which decreases the leading 
edge pressure and pressure gradients. Depending on the flow conditions, as the airfoil 
begins to stall, a vortex forms at the leading edge which grows with time which is 
eventually shed downstream in the wake. This vortex shedding phenomenon alters the 
chordwise pressure distribution on the upper surface of the airfoil resulting in higher 
maximum lift coefficients. Sometimes, however, only a recirculating flow is observed 
at the airfoil trailing edge. As the pitch-up motion progresses the recirculating region 
extends upstream towards the leading edge and eventually a vortex is formed at the 
trailing edge. 
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Rotary wing and fixed wing aircraft designers can enhance the performance of 
the aircraft by taking advantage of the dynamic lift concept. However , the resulting 
undesirable pitching moment variations have deprived helicopters and airplanes of 
the benefits of dynamic lift [Ref. 2]. Carr [Ref. 3] gives a comprehensive review of 
the progress that has been made in the study of the dynamic stall phenomenon. 

Experimental work on steady transonic flow on airfoils was conducted by McDe- 
vitt and Okuno [Ref. 4]. These data cover a wide range of flow conditions and can 
be used for steady code validation. One of the pioneering experimental works on un- 
steady airfoil flows was conducted by McCroskey et al [Ref. 5]. In their experiments, 
unsteady dynamic effects for two-dimensional airfoil flows were studied for the first 
time in conjunction with flows over helicopter blades. They showed the effects of 
unsteady lift and pitching moment on the retreating rotor blades. Their research 
also showed that unsteady motion parameters such as reduced frequency, appeared 
to be more important than airfoil shape in determining the dynamic stall airloads. 
Benchmark experimental data of oscillating and pitching airfoils was also collected by 
Landon [Ref. 6]. The experimental data of references [5] and [6] gives enough quali- 
tative information for unsteady code validation. Recent experiments with oscillating 
airfoils, performed by Chandrasekhara and Brydges [Ref. 7], have shown definitively 
that at subsonic Mach numbers the unsteady flow in the vicinity of the leading edge 
can reach supersonic speeds and generate shocks. 

Along with all the experimental studies being conducted, a great effort is also 
underway to compute unsteady viscous flows. Developments of numerical methods 
for the Navier-Stokes equations [Refs. 8, 9, 10] during the past few years provide 
new tools for the investigation and prediction of airfoil flows. In this investigation 
the compressible thin layer Navier-Stokes equations are solved using a zonal grid 
approach. The objective is to develop a computationally efficient method to study 



two-dimensional unsteady flows. This will be accomplished by developing a solution 
procedure for two grids, an inner viscous grid around the solid body, and an outer 
grid which is coarser representing the outer flow field. The inner grid is free to rotate 
to any angle of attack, and the outer grid remains stationary. 

The idea of moving meshes and dynamic meshes is not new. The dynamic 
and adaptive grid solution method refers to computational grids which are coupled 
to the physical problem that is being solved. An adaptive solution procedure with 
grid points that continually move during the solution process in order to resolve the 
developing gradients has been shown [Ref. 11]. A dynamic type of mesh has been 
applied by Rumsey and Anderson [Ref. 12] to simulate aileron buzz using the thin 
layer Navier-Stokes equations. 

The idea of overlapped or patched grid schemes was used by Rubbert and Lee 
[Ref. 13] with the limitation that the grid lines at the boundaries were continu- 
ous. Benek et al [Ref. 14] developed the “CHIMERA” approach for two-dimensional 
embedded and overlapping grids. They demonstrated that the grid lines at the over- 
lapped boundaries did not have to be continuous, and that the flow quantities could 
be successfully interpolated. The “CHIMERA” approach was later extended to three 
dimensions by Benek et al [Ref. 15] . Rai [Ref. 16] developed a technique for in- 
dependent zonal grids where the flow variables were interpolated across the zonal 
interface. Chesshire and Henshaw [Ref. 17] developed a methodology for solving the 
steady compressible Navier-Stokes equations using multiple overlapped grids. Their 
methodology was demonstrated by solving the governing equations on a composite 
grid that was comprised of an airfoil grid, a leading-edge flap grid, a trailing-edge 
flap grid and a grid for the outer flowfield. Chyu and Davis [Ref. IS] investigated 
unsteady transonic flow using a moving airfoil grid. They developed stationary com- 
putational grids around the airfoil at its lowest and highest angles of attack, and then 
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interpolated a new grid as the airfoil oscillated to intermediate angles-of-attack. Reu 
and Ying [Ref. 19] developed a composite grid approach to study the flow about 
pitching airfoils in a wind tunnel. Their approach consisted of a structured inner grid 
and an unstructured grid for the outer flowfield. 

Unsteady problems have also been solved by oscillating the incoming flow, as 
opposed to moving the grid. However, this method does not work if a solution is 
sought for two or more objects, such as a wing tail combination or a canard wing in 
relative motion to each other. 

B. PURPOSE 

In this work, unsteady compressible flows are investigated using a numerical 
technique which is applied to zonal meshes. The governing equations are solved 
on multiple computational grids, where one of the grids is free to move in unison 
with the solid boundary and the other grid is fixed. The meshes overlap at the 
zonal interface. The scheme that is developed is similar to the “CHIMERA” scheme 
mentioned previously, but we impose the restriction of a circular shape on the zonal 
interface. 

The new approach developed in this study is more computationally efficient 
compared to previous schemes used to study unsteady flows. The zonal grid approach 
avoids the need to regenerate or interpolate entire grids at every angle of attack. The 
zonal interface is in a known position; therefore the entire flowfield does not need to 
be searched. In addition, since the outer grid remains stationary, the metric terms do 
not need to be recomputed at each time step. Another advantage is that the zonal 
interface is circular; therefore the flow variables are interpolated in the circumferential 
direction only. The zonal grid approach also allows for the application of different 
solution methods to the inner and outer grids. For example, viscous solution methods 
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near the solid body and an inviscid method on the outer mesh can be implemented. 

The primary objective of this investigation is to develop and test the zonal grid 
methodology. The space discretization is based on Osher's [Ref. 20] upwind method. 
An implicit scheme is used for time integration. An advantage of upwind schemes is 
that they are naturally dissipative and no explicit artificial dissipation is required. 

In order to meet our goal, a procedure for generating zonal grids was devel- 
oped, along with the zonal flow solver. The effect of grid resolution on the accuracy 
of solutions was studied first for steady flows and the solutions were compared to 
experimental data collected by Harris [Ref. 21]. The unsteady results were validated 
using the experiments by Landon. The dependence of the solution on the location of 
the overlapped zonal region was examined. Viscous steady state solutions were com- 
puted and compared to experimental data. Finally, unsteady flows were investigated 
by computing solutions for airfoils in ramp and oscillatory motion. The ramp motion 
was started at zero angle of incidence, and ramped up to 30 degrees. The computed 
pressure coefficients were verified by comparing with available experimental pressure 
distributions. The boundary layers computed for this case were also compared to an 
interactive boundary layer code [Ref. 22]. The oscillatory test case was verified by 
comparing the computed pressure coefficients with experimental data. This approach 
and the computed results are described in the following chapters. 
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II. GOVERNING EQUATIONS 

In order to compute compressible viscous fluid flow around a body, the conti- 
nuity, momentum and energy equations must be solved simultaneously. The vector 
and the conservation-law form of the compressible Reynolds-averaged Navier-Stokes 
equation is presented. A detailed derivation can be found in [Ref. 23] 

A. CONTINUITY EQUATION 

The continuity equation expresses the conservation of mass law applied to a 
fluid passing through a control volume fixed in space 

^ + (V.pV)=0 (2.1) 

here p is the fluid density and V is the fluid velocity. Equation 2.1 states that the 
net mass flux through a control volume bounding surface must be equal to the time 
rate of change of the mass inside the control volume. In a two-dimensional Cartesian 
coordinate system this equation reads 

dp d d . 

di + Tx (pu) + Tz {pw) = 0 (l2) 

where u and w are velocity components along the x and z directions, respectively. 

B. MOMENTUM EQUATION 

The momentum equation expresses Newton’s second law as applied to a fluid 
element passing through a control volume fixed in space. The momentum equation 
is: 

T~(pV) + V • pVV = pi + V • Ilij (2.3) 
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The first term in equation 2.3 represents the time rate of change of momentum 
per unit volume in the control volume. The second term represents the moment flux 
through the bounding surface of the control volume. / is the body force per unit 
volume and n u is the stress tensor given by 



1 1 ij + P 

where i,j,k = 1,2,3 and 8 i: is the Kronecker delta. 

By substituting equation 2.4 into equation 2.3 and expanding equation 2.3 for 
a two-dimensional Cartesian coordinate system obtain 



du, _ 2 duk 
dx : 3 >J dxk 



12.4' 



0 Qu _ r djz , d_ [2 (c,du _ SuA] , _9_ [ ( dw , 9u \ ] 

Pot PI* dx ^ dx [3 P \~dx dz ) \ + dz [r V dr + dz ) J 



(2.5) 



n Dw _ r dp 1 3 [2 ( cy dw du\] ! 3 ( 3w , du\ 

P~Dt ~ PJ:~ dz + d-z [5/ £ a7 ~ —x)\ + ~x [P (,97 + ~z)\ 

Liquations 2.5 are known as the Navier-Stokes equations for two-dimensional flow. 



C. ENERGY EQUATION 

The energy equation is derived by applying the first law of thermodynamics ( rate 
of change of energy = net heat flux into particle + rate of work done on particle. ) 

-rjj - -p(f x u + fzw) + — (eu + pu-UT xx -wT x: + Q x ) 

+ — (eiC + P W — WT z: — UT i: + Qz) =0 



where e is the total energy per unit volume and Q is the heat addition per unit 
volume. 

The above equations can be rewritten in non-dimensionalized vector form as 



OQ OF dG _ fd Fv dG v \ 
dt ^ dx ^ dz Re \ dx dz J 

where 



( 2 . 6 ) 
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puw 
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. (e + p)u _ 
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T X x = ^P -W z 



Txz — P {u z 4“ U? r ) 

_ 4 / 1 
T ZZ r^P cy^ x 



/a = UT XX + WT XZ + - ; ■■■ —a? 



Pr(7-1) 



<74 = ur x2 + wt zz + 



Pr (7 - 1) 



Re = 



UooL 



where Uoo is the free stream speed and L is a reference length The pressure is related 
to i i, w, p, and e by 

P = { 7 - l)[e - \p{u 2 + u> 2 )] 



In the above equations, the ratio of the specific heats, 7 , is set equal to 1.4 and 
a 2 = ~jp/p is the local speed of sound. 

The density is non-dimensionalized by the free stream density, /><*>, the velocities 
arc non dimensionalized by the free stream speed of sound, and the total energy is 
non dimensionalized by Poo a lo- 
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D. TURBULENCE MODEL 



The Xavier-Stokes equations can completely model fluid flow, however, in or- 
der to resolve turbulent scales at high Reynolds numbers and realistic geometries, 
very high grid densities are required. Therefore, present state-of-the-art algorithms 
and computer technology allows direct simulation of turbulent flows for only simple 
geometries and low Reynolds number flows which are of limited practical interest. 
In order to enable computation of turbulent flows for configurations, of practical in- 
terest, turbulence modeling is used. Turbulence models are implemented with the 
time-averaged forms of the Xavier-Stokes equations. 

Two widely used averaging procedures are the standard time averaging proce- 
dure for incompressible flow, and the mass-averaged approach for compressible flows 
[Ref. 24]. The time averaging procedure destroys high frequency information of the 
turbulence, but the unsteady mean flow information is preserved. 

For the incompressible case, the randomly changing flow variables are replaced 
with their averages plus their fluctuations. For example, in the Cartesian coordinate 
system, the u velocity component is represented as u = u + u', where u is the mean 
velocity and u' is the fluctuation about the mean. The governing equations are time 
averaged, and the average of the fluctuation terms is set equal to zero. 

In the compressible flow case, the mass-weighted variable of the Favre averaging 
approach is used. In this case, u is represented as u = u + u" where u is 

pu 

P 

Here the time average of the doubly primed fluctuating quantities is not equal to zero. 
After the substitution is carried out for all of the fluctuating flow variables, the entire 
equation is time averaged. Next, all the time averaged terms that are doubly primed 
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and multiplied by density are defined to be zero. For example 



The equations of mean motion resulting from the time averaging procedure have 
more unknowns than equations. This constitutes the closure problem of turbulence. 
In order to close these equations a turbulence model is used. Several models have 
been proposed, in this study. 

The Baldwin-Lomax (D — L) turbulence model [Ref. 25] was used. This model 
is a two-layer eddy viscosity model which simulates the effect of turbulence in terms 
of the eddy viscosity coefficient /q. The term y in the stress terms is replaced by 
/< + /q , and the fi/ P r in the heat flux terms is replaced with fij P T + \i t j P Tt . The D — L 
turbulence model is similar to the Cebeci-Smith turbulence model [Ref. 26]. but it 
bypasses the need for finding the edge of the boundary layer by using vorticity instead 
of the boundary layer thickness. This model is adequate for flows which have mild 
pressure gradients, but it is not very suitable for highly separated flows. A complete 
description of the model is given in reference [25]. The basic equations of the model 



In the inner layer, the eddy viscosity is assumed to be proportional to the mixing 
length and vorticity, and in the outer layer it uses an exponentially decaying formula. 



follow. 



The inner eddy viscosity is computed up to the point where it is equal to the outer 
eddy viscosity as shown below. 




(2.7) 



where y is the normal distance from the wall and ^crossover taken at its minimum value 
where it equals y. The inner eddy viscosity is given by: 



(/^)inner — [^ 



inner 
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where 



/ = 



kij 



1 -exp ( -_ 4+) 



'dw_du 



_i_ p w it T y 

y = 

f-^w 



/1 + is an experimentally determined damping constant, k is the Von Karman constant 
The outer eddy viscosity is given by 

/^outer — ^ ^c P PE(i/)\vaI<E-F(j/)klEB- 

^(i/)kleb is the Klebanoff intermittencv factor given by 



F{y )kleb = 
k and C cp are constants. For boundary layers 



l + 5 VbSHi«Y 

V V max / 



F(y)wake — V max^max- 



For wakes and separated boundary layers 



F'(?/)wake — Cxu^y 



max 



U 



2 

DIF 



iv 



The quantity ?/ max is the value of y determined for the maximum value of F max and 
F max is determined by 

F{y) = 2/k’ 

E. TRANSFORMATION TO GENERALIZED COORDINATES 

In order to use an unweighted differencing scheme that facilitate the numerical 
implementation on body fitted coordinate systems suitable for complex configura- 
tions, the equations are transformed to a generalized coordinate system using the 
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following transformations 



c 

s 



x, r 



C = C(i-- 



(2.S) 



The above equations transform the governing equations from the physical domain to 
a body fitted coordinate system. The transformation is carried out by using the chain 
rule of partial differentiation 



d__ f d_ . d_ d_ d_ . d_ 

dx ~^d£ + U d ( ; dz ~ 6 + U dQ 



(2.9) 



For example, the continuity equation would be transformed in the following manner. 

dpu _ c dpu dpv dpv _ c dpu dpv 

dx ^ d{ ^ x d(' dz sc di dc } 

G, s z iCr and G are known as the metrics. The metric terms can determined in the 

following manner. First write the differential expressions for £ and £ 



ds = s xdx + £ z dz; d( = ( X dx + ( Z dz 



( 2 . 11 ) 



In matrix form the expressions become 



■ d£ ' 




' G 6 ■ 


dx 


. d < . 




Cr (z 


dz 



Next we write the differential expression for x and z 



dx = + x^d( 

dz = z^d^ + z<;d( 

In matrix form they become 



dx 




x t 


x < 




dz 




. z t 


z < 


. <*C. 



By comparing equations 2.12 and 2.14 we can write 



1 

H 

yvv 

H 

./■n 


-1 


‘ G G ' 


1 

N 

w 

N* 

1 




. G G . 



(2.13) 



(2.14) 



(2.15) 
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and solving for the metrics we get 



£x = J ~c 

S* = ~JX<1 

C-T = —J~£ 

Cz = J=(, 



(2.16) 



where 



x i 

z < z < : 



(2.17) 



The mapping from (j, 2 ) to (£, £) is one to one if the Jacobian, J . is non singular. 



F. THIN LAYER NAVIER STOKES EQUATIONS 

In order to accurately calculate the normal gradients in the boundary layer, it 
is necessary to make the normal grid spacing very fine close to the solid surfaces, n 
the streamwise direction the flow gradients are not large and no fine grid spacing is 
required. As a result the grid cells near the body have a very high aspect ratio. With 
this type of grid, existing gradients in the streamwise direction are not fully resolved. 
In order to facilitate the numerical implementation of the thin layer approximation 
is employed by retaining only the viscous derivatives normal to the body. The thin 
layer Navier-Stokes equations transformed into a generalized curvilinear system is 
vector form are 



where 



dQ OF dG _ 1 dS 

dt <9£ 6>C Re _<9£ 




P 

pu 

pw 

e 




pU 

puU + £ x p 
pwU + t, x p 
(e + p)U - ZtP . 
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p w 

pu\V + ( x p 
pw\V + C z p 
(e + p)W — ( ( p _ 



The viscous flux term is transformed as 



where 



S 



0 



1 

J 



pmxu^ + (p/3)m 2 C >x 
p.m\W( + (p/3)m 2 ( x 
/imim 3 + (p/3 )m 2 m 4 



m i = £ + C 



m 2 - + CzW C 



m 3 



1^ _d__ ( u 2 +m 2 ) , 1 9g 2 

2 9? 2 i " Pr(-f— 1 ) 9C 



(2- IS) 



m 4 = ( x u + ( 2 w. 

U and W are the contravariant velocity components. These velocity components 
are normal and parallel to the constant £ and ( surfaces, respectively. 



i — st T £x u T 



If" — (t + (x u + C z w 
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III. SOLUTION METHODS 



A. ZONAL GRID GENERATION 

The zonal grids were generated using a software package called GRAPE2D. de- 
veloped at NASA Ames Research Center by Reese L. Sorenson [Ref. 27]. GRAPE2D 
generates field grids by solving the following set of elliptic Poisson equations: 



Six ~ S” = P 
T]xx Vxz Q 



(3.1) 



GRAPE2D was used to generate the inner and outer grid separately. In gener- 
ating the inner grid, the number of circumferential and radial grid points, the spacing 
at the body surface, and the radius and shape of the outer boundary were specified. 
Once the inner grid was generated, the outer boundary of the inner grid was used as 
the inner boundary for the outer grid. 

Due to the placement of the reference axis at the quarter chord of the airfoil, 
the grid spacing in front of the airfoil is larger than aft of the airfoil. This causes a 
problem in specifying an initial grid spacing for the outer grid. In order to obtain 
the smoothest overall transition from the inner grid to the outer grid, the starting 
spacing for the outer grid was obtained by averaging the minimum and maximum 
spacing at the outer boundary of the inner grid. The outer boundary of the outer 
grid was prescribed as a rectangular region, usually about six chord lengths away 
from the body. 

The perfectly circular zonal boundary was accomplished by reading the over- 
lapped regions x and 0 coordinates and finding an average radius which was measured 
from the center of the grid, regardless of the location of the inner body. Next, an 
arbitrary grid point was chosen to use as a reference point. Then an angle theta, 
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6 , was calculated which represented the angle between the reference grid point and 
the grid point under consideration. Then having the average radius and the angle 0 
for the grid points in the overlapped region, a new set of x and y coordinates was 
calculated using the following equations: 

x = r cos 0 , y = r sin 0 (3.2) 

The overlapping of the two grids was performed by adding the next to last layer 
of grid points from the inner grid to the outer grid, and adding the second layer of 
grid points from the outer grid to the outer boundary of the inner grid. This resulted 
in an overlapped region of three grid points or two grid cells in the radial direction. 
The grid was also overlapped by a grid point in the circumferential direction. 



B. NUMERICAL SCHEME 

The numerical integration is performed using an upwind-biased, factorized, 
iterative, implicit numerical scheme [Ref. 20] given by 



I + h ( (V\A* k + A'.4- t ) 

r + MVJBJ + ±{B- t - Re-'6 ( M l , k )Y X (Q'f - Q' t ) = (J 
-I (Ql„ - Qlk) + M^+ w - FL m ,„) 

+'k(Gf. l+l/2 - - Re-'h(Sl um - S? t _„ 2 )] 

In equation 3.3, h $ = Ar/A£, etc., = ( dF/dQ ), etc., are the flux Jacobian 
matrices, and A, V, and 6 are the forward, backward and central difference operators, 
respectively. The quantities Ei+ 1 / 2 ,/:, C/^+ 1 / 2 , and Si } k+ 1/2 are numerical fluxes. The 
inviscid fluxes F and G are evaluated using Osher's upwinding scheme [Ref. 19]. The 
numerical fluxes for a third-order accurate scheme are given by 



C + l/2,fc - C+l/2,fc + | ,k + -^^1 + 1/2. 

6 -^^i+3/2, k + “‘“^1 + 1/2 , k = F{Qi,k, Qi+l,k) + 
l [AF+(Q 1+U , Q itk ) + 2A F+{Q iJe , Q t+hk )) - 
I (AF (Ci.fc, Qi+\,k) + 2AF (C1+1,/:, Ci./c)] — 



(3.4) 
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\\ here F is the first-order numerical flux for the explicit Osher’s scheme given by 



-R 1 / 2 , k 



fU + f,+ f Q,t ' {F,+ - f;} , IQ 

J Qx 



where F 7 = + F~ , F± — (^) . and A F ± are the corrections to obtain a 

higher order accuracy. The Oslier scheme evaluates the flux assuming a shock tube 
solution where F q is piecewise continuous and yields good predictions of the flux, 
especially at supersonic Mach numbers. For the linearization of the left-hand side of 
Eq. 3.3, the flux Jacobian matrices A, B are evaluated by the Steger- Warming [Ref. 
28] flux-vector splitting. 

Time accuracy of the implicit numerical solution is obtained by performing 
Newton iteration to convergence for each time step. The approximation of Q n+l 
at each subiteration is the quantity Q p . When p > 2, during a given subiteration. 
Q p = 0 n+1 i but when p = 1 and no subiteration is performed, then Q p = Q n , and 
Q p+1 = Q n+1 ■ By subiterating to convergence, linearization and factorization errors 
are minimized, because the left-hand side of Eq. 3.3, can be driven to zero at each 
time step. 

The linearization errors are eliminated by subiteration to convergence. Typ- 
ically, two to three subiterations are sufficient to drop the residuals two orders of 
magnitude during the Newton iteration process. 

High order accurate shock-capturing schemes have some limitations; they may 
select a nonphysical solution, they produce spurious oscillations and they may de- 
velop a nonlinear instability in nonsmooth and discontinuous flow regions. More ap- 
propriate high-order shock-capturing schemes, suitable for the computation of weak 
solutions are the TVD schemes, described in detail in [Ref. 29, 30]. In the present 
study, the Osher-Chakravarthy [Ref. 31] TVD scheme is used. This TVD scheme has 
flux limiters which impose constraints on the gradients of the fluxes. The flux-limited 
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values A f ± are computed from the unlimited fluxes A/ ± as follows 

-Vi+3/2,* = rninmocl 
A /j+i/ 2 ,(t = minmod 
A/i+i/ 2 ,jfc = minmod 

A/ 7 - 1/2 t = minmod 

where the minmod operator is defined by 

mimno(/[r. y] = sign(x) x mar[0, mm{|i|, ysfgn(x)}] (3.6) 

The viscous duxes 5 1( t + i/2 are computed with central differences as follows 



A/ i +3/2, fc' ^A/- 1/2 , fc 

A/ i+l/2,/:’ J-^fi+3/2,k 

A fZwMt-i/i.k 
A/i-l/2,t-‘ r ^A/7i/2,fc 



, A: 4-1/2 — S[Q ty k+l/2, (Qc)ufc+l/2> Ci.Jt+1/2] 


(3.7) 


Q x,k + l /2 — 7) {Qi,k “1“ Qi,k — l) 


(3.3) 


(Qc)*\k+l/2 — Qi,k+l Qiyk 


(3.9) 



The experimental Reynolds numbers based on the chord length for the test cases 
examined are in the range /?e c % 3.0 x 10 6 — 5.0 x 10 6 , and it is expected that the how 
is mostly turbulent. Transitional dow is expected to have a small effect at regions very 
close to the leading edge. Present knowledge about boundary layer transition does 
not enable computation and modeling of the transition regime. Therefore, only fully 
turbulent solutions were computed. In the present work, the widely used two-layer 
BaUlwin-Lomax turbulence model was used. The effectiveness of other turbulence 
models, such as the Johnson-I\ing model [Ref. 32] and the RNG based algebraic 
model [Ref. 33] for steady and unsteady dows, was investigated in references [34] and 

[35]. 

C. BOUNDARY CONDITIONS 

The solutions on the two grids are computed separately, with the inner and 
outer solutions communicating through the zonal interface boundary. The inner grid 
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surrounds the airfoil and includes the boundary layer region and the wake for viscous 
solutions. Inviscid solutions are obtained by applying the flow tangency slip condition 
where the normal contravariant velocity component. IV’, is set equal to zero on the 
surface. For viscous solutions, the nonslip condition is applied for the velocities on the 
airfoil surface. In both cases the density and pressure are obtained from the interior 
grid points by simple extrapolation. For unsteady solutions, the surface velocity is 
set equal to the airfoil speed obtained by the prescribed airfoil motion as follows 

^ — W = y (stCx 

Unsteady solutions for pitching and oscillating airfoils are obtained by rotating the 
inner grid only. Therefore, only the metrics of the inner grid must be reevaluated for 
each time step. 

At the inner zonal interface, the flow variables are obtained from the interior of 
the outer grid solution. Similarly, the inner zonal boundary of the outer grid obtains 
boundary information from the interior of the inner grid. The inner and outer grid 
radial lines are not aligned, in general. The relative location of the two grids with 
respect to the inertial reference frame as the inner grid rotates is computed. These 
distances between neighboring points at the zonal interface are used for a weighted 
averaging of the conservative variables. 

All flows were computed for subsonic free-stream speeds. At the subsonic inflow 
and outflow boundaries of the outer grid, the flow variables were reevaluated using 
zero-order Riemann invariant extrapolation. At the inflow boundary, there are three 
incoming and one outgoing characteristics. Therefore, three variables, the density p, 
the normal velocity w , and the pressure p, are specified and the fourth variable, the 
axial velocity u is extrapolated from the interior. The inflow boundary conditions are 
given by 
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p 1 = 



(!/-,-!) 



- S 1 



Py 
P > o 



«i = ^(Rt-Ri) 

a i — (R+ + R 2 )!'l 



1 3. 10) 



tUl = Wy, 

Pi =(*?) 

1 km e /?+, Ro are the incoming and outgoing Riemann invariants given by 



Rt = u + 2a.x>/(7 - !)• ^2 = u 2 - 2a 2 /b - 1) 

\t the outflow boundary there are one incoming and three outgoing characteristic*. 

I lierefore only one quantity, the pressure, is specified, while the others are extrap- 
olated from the interior. For the density and normal velocity, simple first-order 
extrapolation is used, and the axial outflow velocity is obtained from the zero ordei 
outgoing Riemann invariant. The outflow boundary conditions are given by 



P i = P 2 

= Rt -2ai/0 - l) 

a i = \pip\Tpi 



o. m 



ir ! — tt’9 



Pi = /'x 
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were solved. For the unsteady flow solutions, the outer grid remained stationary and 
the metrics were not reevaluated at each time step. 
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IV. RESULTS AND DISCUSSION 



The validity of the zonal grid approach was first investigated for inviscid solu- 
tions. An advantage of the present approach is that different grid densities may be 
used for the inner and outer grids. However, the accuracy and the conservative char- 
acter of the solution for different locations of the zonal interface and grid densities 
must be assessed. 

First, the accuracy of the computed results for different inner and outer grid 
densities was evaluated. The effect of the location of the zonal interface relative to the 
airfoil on the accuracy of the solution was also investigated. Then viscous solutions at 
fixed angles of incidence, up to approximately the static stall angle, were computed. 
Finally, unsteady flow responses to a ramp motion at subsonic free-stream speed of 
iV/oo = 0.3 and for an oscillation at a free-stream speed of M,^, = 0.6 were computed. 

A. STEADY STATE SOLUTIONS 
1. Inviscid Test Cases 

Preliminary test cases were computed using coarse meshes with an inviscid 
solution. A two-block grid consisting of an SI x 40 point O-type inner grid and an 
SI x 22 point O-type outer grid was used as a baseline grid for the inviscid solutions. 
Table 4.1 gives the inviscid grids that were tested. 

a. Case 1. Baseline Grid: 81x40 Inner and 81x22 Outer 

An inviscid solution using the baseline grid for subsonic flow over a 
XACA-0012 airfoil at = O.S, a = —.1° was obtained. The baseline grid is given in 
figure 4.2. The distribution of the computed surface pressure coefficient is compared 
with the measurements of [Ref. 4] in Fig. 4 . 1 . Agreement with the experimental data 
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TABLE 4.1: GRID DENSITIES OF THE INVISCID SOLUTIONS COMPUTED 
FOR A NACA 0012 AIRFOIL AT = .S AND a = -0.1. 



Case 


Inner 


Grid 


Outer 


Grid 


Inner Grid Radius 


i 


SI 


X 


40 


SI 


X 


22 


1.5 


x chord 


2 


SI 


X 


40 


41 


X 


22 


1.5 


x chord 


3 


SI 


X 


40 


SI 


X 


22 


1.0 


x chord 


4 


SI 


X 


40 


SI 


X 


22 


.75 


x chord 


5 


SI 


X 


IS 


SI 


X 


19 


Oi 


mlGrid 


6 


SI 


X 


20 


41 


X 


12 


1.5 


x chord 



is satisfactory for an Euler solution. It can be seen that with this mesh the strength 
of the shock is captured, but the location is lagged by 5 percent of the chord. The 
baseline solution is converged at 2000 iterations. 



a = -0.1 M = 0.8 




Figure 4.1: M ^ = 0.8 Computed Solution With Baseline Grid Compared to Exper- 
iment. 
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Figure 4.2: Inner Grid 81 x 40 Outer Grid is 81 x 22 
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b. Case 2. Baseline Inner Grid With a 41x22 Outer Grid 



Next the effect of the outer grid density is investigated. The grid for 
Case 2 was generated by starting with the baseline grid and removing every other grid 
point from the outer grid. The resulting grid is given in Fig. 4.4. The Euler solution 
is compared to experiment in Fig. 4.3. Overall the converged solution for this case 
agrees with the baseline solution except it is noticed that the shock location for the 
upper and lower surface are slightly farther apart than the baseline grid solution. 



a = -0.1 M = 0.8 




Figure 4.3: = 0.S Computed pressure coefficient solution compared to experi- 

mental data. 



25 




Figure 4.4: Inner Grid 81 x 40 Outer Grid is 41 x 22 
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c. Case 3. Overlap Boundary Set to 1 Chord Length Away From 
Body 

In this case the effect of the overlap boundary location is studied. The 
grid generated is given in Fig. 4.6. It is seen that the transition from the inner to 
the outer grid is not as smooth as for the cases where the boundary was located 50 
percent farther from the airfoil. At the leading and trailing edges, the overlapped 
boundary is only a half chord length away. In figure 4.5 the converged solution is 
given. With this grid very little effect on the shock strength and location is seen. The 
pressure coefficient before the shock is slightly underpredicted and after the shock is 
slightly overpredicted. 



a = -0.1 M = 0.8 




Figure 4.5: M ^ = 0.8 Computed Solution Compared to Experiment. 
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Figure 4.6: Inner Grid 81 x 40 Outer Grid is 81 x 22 
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d. Case 4. Overlap Boundary Set to .75 Chord Length Away 
From Airfoil 

In order to further investigate the tendencies observed in Case 3, an- 
other grid is developed with the zonal interface closer to the airfoil. The overlapped 
region is only a quarter chord away from the leading and trailing edges of the airfoil 
as seen in Fig. 4.S. The solution obtained is compared to experiment in Fig. 4.7. 
The pressure coefficient is again slightly underpredicted before the shock and slightly 
overpreclicted after the shock. It is also observed that the zonal interface location 
relative to the airfoil has little effect on the solution. 



a = -0.1 M = 0.8 




Figure 4.7: M ^ = 0.8 Computed Solution Compared to Experiment. 
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e. Case 5. Grid With Oval Interface 



I he solution obtained for cases 3 and 4 showed that the location of 
the zonal interface had very little effect on the computed solutions. The ability of 
the zonal interface to pass flow discontinuities was also studied. It was known from 
the previous cases that the shock location was near the the mid chord point. The 
inner grid for this test case had to be generated so that the overlapped region was 
very close to the upper and lower surface of the airfoil. This was accomplished by 
generating an oval inner grid and an oval zonal interface. This grid is shown in figure 
4.10. The computed surface pressure coefficient distribution is shown in figure 4.9. It 
is in agreement with the experimental data and with the previous computed solutions. 
The computed flow quantities, such as density and pressure, showed that the zonal 
approach used can pass shocks through the zonal interface. Figure 4.11 shows Mach 
contours which smoothly pass through the overlapped boundary. 



a = -0.1 M = 0.8 




Figure 4.9: M ^ = 0.8 Computed Solution Compared to Experiment. 
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Figure 4.10: Inner Grid SI x 18 Outer Grid is SI x 19 
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Figure 4.11: Mach Contours: Overlapped Region's Ability to Pass Flow Disconti- 
nuities. 
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f. Case 6. Baseline Grid With Half the Radial Grid Points On 



Outer Grid 

In the previous cases, the effect of circumferential resolution was stud- 
ied. Next, the effect of the radial resolution on the computed solutions is investigated. 
Figure 4.13 shows the grid generated to test these effects. The computed surface 
pressure coefficient distribution (in figure 4.12) was comparable with the solutions 
obtained with denser outer grids. Computations with an even coarser grid, e.g., a 
41 x 11 point grid, predicted the shock location even further downstream due to lack 
of streamwise resolution. 



a = -0.1 M = 0.8 




Figure 4.12: = 0.8 Computed Solution Compared to Experiment. 
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2. Viscous Test Case 



The Euler solutions presented in the previous section were used to study 
the effects of different grid densities on the computed solutions. Having confidence 
in the predictions of the code, the next step was to generate a viscous grid with the 
airfoil's quarter chord point at the center of the inner grid. 

Viscous, subsonic flow solutions were obtained at the following fixed angles 
of attack: 3.27°. 4.97°, 6.69°, 8.38°, 9.27°, 10.12°, 10.99° and 11.90°. The flow condi- 
tions of the measurements reported in [Ref. 21] were used, e.g., 3/^ = 0.3. Re = 
4.0 x 10 6 . The spacing of the grid was set to 0.0005 for the first grid point above 
the surface. The quarter chord point of the airfoil was set at the center of rotation 
so the airfoil could be ramped about the quarter chord point. These solutions were 
obtained on a 181 x 56 point viscous inner grid and a 181 x 26 point inviscid outer 
grid shown in figure 4.14. 

Solutions were also computed on a grid with half the streamwise resolution, 
e.g.. a 91 x56 point grid. The computed surface pressure coefficient distributions using 
the 181 x 56 inner grid and 181 x 26 outer grid, are compared to experimental results 
in figures 4.15 through 4.23. 

Solutions for fixed angles of incidence were obtained by two methods. First 
by rotating the inner grid to the specified angle of incidence and setting the oncoming 
flow to zero degrees. Second by rotating the flow to the angle of incidence and 
leaving the inner grid at zero angle relative to the outer grid. The computed pressure 
coefficients and boundary layers were the same for both cases. 

For the low Mach number viscous solutions, no flux limiting was applied. 
It is seen that the computed results closely agree with the experimental data. At the 
higher angles of incidence the suction peak is not exactly captured. This is probably 
due to lack of grid resolution at the leading edge of the airfoil. 
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Figure 4.14: A NACA-0012 Airfoil Grid Centered at the Quarter Chord Point. 
IS 1x58 Inner Grid and 181x26 Outer. 
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a=3.27 deg., M=0.30, Re=2700000 




Axial Location x/c 



Figure 4.15: Viscous Computed Solution Compared to Experiment for a = 3.27°. 



a=4.97 deg., M=0.30, Re=2700000 




Figure 4.16: Viscous Computed Solution Compared to Experiment for a = 4.97°. 
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a=6.69 deg., M=0.30, Re=2700000 




Figure 4.17: Viscous Computed Solution Compared to Experiment for a = 6.69° 



a=7.54 deg., M=0.30, Re=2700000 




Axial Location x/c 



Figure 4.18: Viscous Computed Solution Compared to Experiment for a = 7.54° 
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a=8.38 deg., M=0.30, Re=2700000 




Figure 4.19: Viscous Computed Solution Compared to Experiment for a = 8.38°. 



a=9.27 deg., M=0.30, Re=2700000 




Axial Location x/c 



Figure 4.20: Viscous Computed Solution Compared to Experiment for a = 9.27°. 
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a=l0.12 deg., M=0.30, Re=2700000 
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Figure 4.21: Viscous Computed Solution Compared to Experiment for a = 10.12 
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a=10.99 deg., M=0.30, Re=2700000 




Axial Location x/c 



Figure 4.22: Viscous Computed Solution Compared to Experiment for a = 10. 99°. 
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Figure 4.23: Viscous Computed Solution Compared to Experiment for a = 11.99°. 
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B. RAMP MOTION SOLUTION 



The unsteady solution for a rainp motion from a = 0° to q = 30° at = 
0.3, Re = 2.7 x 10 e and reduced frequency , k = 0.0127 was obtained on both a 
91 x 56 point inner grid and a 181 x 56 point inner grid. The pitch rate for the ramp 
motion k is defined as k = ac/'lb '^ . The computed lift response is compared with 
the experimental measurements of [Ref. 6] in Figure 4.24. 



a(t)=0.0 to 15.5 deg., M=0.30, k=0.0l27, Re=2700000 




Figure 4.24: Comparison of the Measured and Computed Lift for the Ramp Motion 

Both the coarse and fine grid solutions closely predict the measured lift. How- 
ever, at the higher angles of attack, the finer grid gives higher lift. The computed 
surface pressure coefficient distributions at several angles of incidence are compared in 
Figures 4.25 - 4.38. Experimental surface pressure coefficients were available for the 
following angles of incidence: 2.94°, 5.84°, 8.91°, 11.76°, 15.5°, and they are displayed 
as diamonds in the figures. The computed surface pressure coefficient distribution 
is in good agreement with the measured data over the entire incidence range. The 
computed flowfield at the maximum angle of incidence, a = 15.5°, is mostly attached. 
A small separated flow region exists at the trailing e .ge region only. At a higher angle 
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of attack, a % 17.0°, the computed solution shows the development of the dynamic 
stall vortex in the leading edge region. 
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Figure 4.25: Computed ramp solution at a = 1.S6 0 . 
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Figure 4.26: Computed ramp solution compared to experimental data at a = 2.94°. 
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a=4.14 M=0.3 k=.0127 Re=2700000 




Figure 4.27: Computed ramp solution at a = 4.14° 



a=4.87 M=0.3 k=.0127 Re=2700000 




Figure 4.28: Computed ramp solution at a = 4.87° 
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a=5.85 M=0.3 k=.0127 Re=2700000 




Figure 4.29: Computed ramp solution compared to experimental data at a = 5.85°. 



a=6.72 M=0.3 k=,0127 Re=2700000 




Figure 4.30: Computed ramp solution at a = 6.72°. 
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a=8.02 M=0.3 k=.0127 Re=2700000 




Figure 4.31: Computed ramp solution at a = 8.02°. 



a=8.91 M=0.3 k=.0127 Re=2700000 




Figure 4.32: Computed ramp solution compared to experimental data at a = 8.91 \ 
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a=9.85 M=0.3 k=.01 27 Re=2700000 




Figure 4.33: Computed ramp solution at a = 9.85°. 



a=10.80 M=0.3 k=.01 27 Re=2700000 




Figure 4.34: Computed ramp solution at a = 10.80°. 
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a=1 1.77 M=0.3 k=.0127 Re=2700000 




Figure 4.35: Computed ramp solution compared to experimental data at 

a = 11.77°. 



a=12.84 M=0.3 k=.0127 Re=2700000 




Figure 4.36: Computed ramp solution at a — 12.84°. 
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a=13.89 M=0.3 k=.0127 Re=2700000 




Figure 4.37: Computed ramp solution at a — 13.S9 0 . 



a=l5.55 M=0.3 k=.0127 Re=2700000 




Figure 4.38: Computed ramp solution compared to experimental data at 

a = 15.55°. 
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1. Boundary Layer Comparisons For Ramp Motion 



Early on in the development of the software, the computed pressure coef- 
ficients agreed quite well with the experimental data. As test cases at higher angles- 
of-attack were tested, it was discovered that the flow was not separating from the 
trailing edge of the airfoil. The problem was discovered by invest igating the boundary 
layer profiles. A problem in the turbulence model was discovered and easily fixed. 
The lesson here is that although the computed and experimental pressure coefficients 
agree quite well it is important to check that the boundary layer profiles are reason- 
able. In figures 4.39 through 4.48 the computed boundary layers are compared to 
boundary layers computed by the interactive viscous inviscid boundary layer method 
of reference [22]. 

The comparisons are made at x/c = .5 and at x/c = .9 for angles-of- 
attack ranging from 2.94° to 15.5°. The computed boundary layer profiles compare 
quite well up to 15.5°. At x/c = .9 for the angle-of-attack of 15.5°, the comparisons 
diverge. This happens because of the different turbulence models used. At this angle- 
of attack, the flow is separating at the trailing edge so the turbulence models used 
become very important. 
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a=2.94 x/c=.5 M=0.3 k=.0127 Re=2700000 




Figure 4.39: Comparison of computed boundary layer with an interactive boundary 
layer program at the 90% chord for a = 2.94 degrees. 



a=2.94 x/c=.9 M=0.3 k=.0127 Re=2700000 




Figure 4.40: Comparison of computed boundary layer with an interactive boundary 
layer program at the 90% chord for a — 2.94 degrees. 



52 



a=5.85 x/c=.5 M=0.3 k=.0127 Re=2700000 




0 0.2 0.4 0.6 0.8 1 

Normalized Velocity u/U 



Figure 4.41: Comparison of computed boundary layer with an interactive boundary 
layer program at the 50% chord for a = 5.85 degrees. 




0 0.2 0.4 0.6 0.8 1 

Normalized Velocity u/U 

Figure 4.42: Comparison of computed boundary layer with an interactive boundary 
layer program at the 90% chord for a = 5.85 degrees. 
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a=8.91 x/c=.5 M=0.3 k=.0127 Re=2700000 




Figure 4.43: Comparison of computed boundary layer with an interactive boundary 
layer program at the 50% chord for a = 8.91 degrees. 




0 0.2 0.4 0.6 0.8 1 

Normalized Velocity u/U 

Figure 4.44: Comparison of computed boundary layer with an interactive boundary 
layer program at the 90% chord for a = 8.91 degrees. 
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0 0.2 0.4 0.6 0.8 1 

Normalized Velocity u/U 



Figure 4.45: Comparison of computed boundary layer with an interactive boundary 
layer program at the 50% chord for a = 11.77 degrees. 




0 0.2 0.4 0.6 0.8 1 

Normalized Velocity u/U 



Figure 4.46: Comparison of computed boundary layer with an interactive boundary 
layer program at the 90% chord for a = 11.77 degrees. 
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a=15.55 x/c=.5 M=0.3 k=.0127 Re=2700000 




Figure 4.47: Comparison of computed boundary layer with an interactive boundary 
layer program at the 50% chord for a = 15.55 degrees. 



a=15.55 x/c=.9 M=0.3 k=.0127 Re=2700000 




Figure 4.48: Comparison of computed boundary layer with an interactive boundary 
layer program at the 90% chord for a — 15.55 degrees. 
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2. Computed Ramp Flow Details 



In this section, flow features of the computed ramp motion are shown 
I liese features are demonstrated by the density contour lines. Mach contour lines, 
mticity contour lines and mass-flux contours. The contour lines represent areas of 
i lie flowfield. for which the parameter of interest, remains constant. Clustered contour 
lines represent areas where the parameter is rapidly changing. 

In figures, 4.49 through 4.64. density. Mach number, vortieity and ma» 
Mux contours are displayed for a = 1.86° to a = 17.60°'. In figures 4.65 through 4.1 II) 
i lie Mach number contours are replaced by pressure contours. This was done beeau-e 
i In- strong gradients, caused by the vortices, tend to not give clear Mach number 
information. 

In figure 4.63. a = 16.50°. the beginning of a vortex is visible near Un- 
loading edge. At a = 17.00°, in figure 4.64 the vortex has moved 5/7 of chord 
downstream, and another vortex is starting at the leading edge. At o = 19. DO 
the first two vortices merge into one stronger vortex located near the 50 % chord 
point. Two smaller vortices are also clearly visible in the mass-flux contours of liiiure 
1.08. The original vortex now detaches from the surface of the airfoil at 19.50 . At 
n = 20.50°, a second vortex is shed from the leading edge area, and a the vortex shown 
at 19.5° has moved to the trailing edge. At a — 22.50°. the trailing edge vortex ha-* 
been shed downstream. The whole process seems to repeat itself at a faster rate, , 1 - 
can Ire seen in the remaining figures as the airfoil ramps up to 30°. The dynamic -tall 
ulienomenon is observed in terms of density, pressure and vortieity. The nose-down 
pitching moment is caused by the vortex sitting at the trailing edge. It is important 
i In- i lie vortices do not dissipate and do not get distorted as they pass through the 
/oual interface. 



Next the solution was continued to angles beyond stall. Traces of alternate 
vortex shedding from the trailing edge can be seen at a % 27°. At high angles of 
incidence the alternating vortex shedding is very well demonstrated in figures 4.91 
through 1.110. These computed solutions are in general agreement with the findinu» 
of tin? experimental investigations of C’hanclrasekahara et al. Also it can be seen that 
the zonal method developed is capable of computing unsteady flows at very high 
uugles-of-attack showing at least qualitative agreement with the experiment. 
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DENSITY MACH NUMBER 







Figure 4.49: Ramp Motion Flow Details, = .3, k = .0127, Re = 2.7 x 10 b . 
a = 1.86°. 
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Figure 4.50: Ramp Motion Flow Details, M 0 0 = .3, k = .0127, Re = 2.7 x 10 b . 
a = 2.94°. 
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Figure 4.51: 
a = 4.14°. 



Ramp Motion Flow Details, M oo 



= .3, k = .0127, Re = 2.7 x 10 6 , 
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Figure 4.52: Ramp Motion Flow Details, = .3, k = .0127, Re = 2.7 x 10 6 , 
a = 5.85°. 
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Figure 4.53: Ramp Motion Flow Details, A/co = -3, k — .0127, Re = 2.7 x 10 6 . 
a = 6.72°. 
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Figure 4.54: Ramp Motion Flow Details, = .3, k = .0127, Re = 2.7 x 10 b . 
a = 8.02°. 
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Ramp Motion Flow Details, M 0 0 = .3, k = .0127, Re = 2.7 x 10 6 . 



Figure 4.55: 

a = 8.91°. 
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Figure 4.56: Ramp Motion Flow Details, M&, = .3, k = .0127, Re = 2.7 x 1 0 6 , 
a = 9.85°. 
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Figure 4.57: Ramp Motion Flow Details, = .3, k = .0127, Re = 2.7 x 10°, 
a = 10.80°. 
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Ramp Motion Flow Details, = .3, h = .0127, Re = 2.7 x 10 6 , 



Figure 4.58: 
a = 11.77°. 



68 



DENSITY MACH NUMBER 







Figure 4.59: Ramp Motion Flow Details, M 0 0 = .3, k = .0127, Re = 2.7 x 10 6 . 
a = 12.84°. 
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Ramp Motion Flow Details, M 0 0 = .3, k = .0127, Re = 2.7 x 10 6 , 



Figure 4.60: 

a = 13.39°. 
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Ramp Motion Flow Details, M ^ = .3, k = .0127, Re — 2.7 x 10 tj . 



Figure 4.61: 

a = 15.55°. 
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Figure 4.62: Ramp Motion Flow Details, M ^ = .3, k = .0127, Re = 2.7 x 10 b , 
a = 16.00°. 
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Figure 4.63: Ramp Motion Flow Details, M 0 0 = .3, k = .0127, Re — 2.7 x 1 0 e , 
a = 16.50°. 
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Figure 4.64: Ramp Motion Flow Details, M oo = .3, k — .0127, Re = 2.7 x 10 6 , 
Q = 17.00°. 
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Figure 4.65: Ramp Motion Flow Details, Moo = -3, k = .0127, Re = 2.7 x 10 6 . 
a = 17.50°. 
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Ramp Motion Flow Details, M ^ = .3, k = .0127, Re = 2.7 x 10 6 , 



Figure 4.66: 
q = 18.00°. 
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Figure 4.67: Ramp Motion Flow Details. = .3, k = .0127, Re = 2.7 x 10 6 . 
a = 18.50°. 
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Figure 4.68: Ramp Motion Flow Details, M 0 0 = .3, k = .0127, Re = 2.7 x 10 c , 
q = 19.00°. 
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Figure 4.69: Ramp Motion Flow Details, = .3, k = .0127, Re = 2.7 x 10 b . 
a = 19.50°. 
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Figure 4.70: Ramp Motion Flow Details, = .3, k — .0127, Re = 2.7 x 10°. 
a = 20.00°. 
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Figure 4.71: Ramp Motion Flow Details, M ^ = .3, k = .0127, Re = 2.7 x 10 6 . 
q = 20.50°. 
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Figure 4.72: Ramp Motion Flow Details, A/^ = .3, k = .0127, Re = 2.7 x 10 e , 
q = 21.00°. 
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Figure 4.73: Ramp Motion Flow Details, = .3, k = .0127, Re = 2.7 x 1 0 6 , 
a = 21.50°. 
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Figure 4.74: Ramp Motion Flow Details, M ^ — .3, k = .0127, Re — 2.7 x 10 u . 
q = 22.00°. 
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Figure 4.75: Ramp Motion Flow Details, Moo = -3, k = .0127, Re = 2.7 x 10°, 
a = 22.50°. 
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Figure 4.76: Ramp Motion Flow Details, = .3, k = .0127, Re = 2.7 x 10 b , 
a = 23.00°. 
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Figure 4.77: Ramp Motion Flow Details, M oo = .3, k = .0127, Re = 2.7 x 1 0 6 , 
q = 23.50°. 
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Figure 4.78: Ramp Motion Flow Details, M <*, = .3, k = .0127, i?e = 2.7 x 10 G , 
a = 24.00°. 
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Figure 4.79: Ramp Motion Flow Details, = .3, k = .0127, Re = 2.7 x 10 b , 
q = 24.50°. 
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Figure 4.80: Ramp Motion Flow Details, — .3, k = .0127, Re = 2.7 x 10°, 
a = 25.00°. 
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Figure 4.81: Ramp Motion Flow Details, = .3, k = .0127, Re = 2.7 x 10 b , 
a = 25.50°. 
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Figure 4.82: Ramp Motion Flow Details, M^ = .3, k = .0127, Re = 2.7 x 10". 

a = 26.00°. 
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Figure 4.83: Ramp Motion Flow Details. = .3, k — .0127. Re — 2.7 x 10 b . 
a = 26.50°. 
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Figure 4.84: Ramp Motion Flow Details, M oo = .3, k = .0127, Re — 2.7 x lO 6 . 
o = 27.00°. 
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Figure 4.85: Ramp Motion Flow Details, = .3, k = .0127, Re = 2.7 x 10°. 
a = 27.50°. 
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Figure 4.86: Ramp Motion Flow Details, M*, = .3, k = .0127, Re = 2.7 x 10°, 
a = 28.00°. 
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Figure 4.87: Ramp Motion Flow Details. M ^ = .3, k — .0127. Re = 2.7 x 10 b . 
n = 28.50*. 
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Figure 4.88: Ramp Motion Flow Details. — .3, k = .0127, Re = 2.7 x 10 b . 

n = 20.00°. 
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Figure 4.89: Ramp Motion Flow Details, = .3, k = .0127, Re = 2.7 x 10°. 
n — 29.50°. 
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Figure 4.90: Ramp Motion Flow Details. M iyz , = .3, k = .0127, Rt = 2.7 x I0 b . 

n = 20.00°'. 
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Figure 4.91: Ramp Motion Flow Details, — .3, k — .0127, Re = 2.7 x I0 b . 
n = 31.00°. 
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Figure 4.92: Ramp Motion Flow Details. = .3, k — .0127, Re = 2.7 x 10 G . 
n = 32.00°. 
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Figure 4.93: Ramp Motion Flow Details. ,\L X = .3, k = .0127. Re = 2.7 x 10 b . 
n = 33.00°. 
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Figure 4.94: Ramp Motion Flow Details, = .3, k = .0127, Re = 2.7 x 10°, 
a = 34.00°. 
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Figure 4.95: Ramp Motion Flow Details, 4/^ = .3, k = .0127, Re = 2.7 x 10 b . 
n = 35.00°. 
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Figure 4.96: Ramp Motion Flow Details, = .3, k = .0127, Re = 2.7 x 10 6 , 
o = 36.00°. 
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Figure 4.97: Ramp Motion Flow Details. M ^ = .3, k = .0127. Re = 2.7 x 10 b . 
«. = 37.00°. 
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Figure 4.98: Ramp Motion Flow Details, = .3, k = .0127, Re = 2.7 x 10 6 , 
a = 38.00°. 
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Figure 4.99: Ramp Motion Flow Details. ,\[ 00 = .3, k = .0127, Re = 2.7 x 10 b . 
.. = 39.00°. 
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Figure 4.100: Ramp Motion Flow Details. M <*, = .3, k = .0127, Re = 2.7 x 10F 

n = 10.00°. 
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Figure 4.101: Ramp Motion Flow Details, A/cc = .3, k = .0127, Re = 2.7 x 10 6 . 
a = 41.00°. 
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I igure 4.102: Ramp Motion Flow Details. .1/^ = .3, k = .0127, Re = 2.7 x 10F 
•. = 12 . 00 °. 
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Figure 4.103: Ramp Motion Flow Details. M ^ = .3, k = .0127, Re = 2.7 x 10*'. 
n = 43.00°. 
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Figure 4.104: Ramp Motion Flow Details. = .3 , k = .0127, Re = 2.7 x 10'’. 

< i = 14.00*. 
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DENSITY PRESSURE 






Figure 4.105: Ramp Motion Flow Details. ,\[ 00 = .3, k = .0127. Re = 2.7 x 10''. 
m = lo.00°. 



OENSITY PRESSURE 




Figure 4 . 106 : Ramp Motion Flow Details, = .3, k = .0127, Re = 2.7 x 1 0 6 , 
a = 46.00°. 
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DENSITY PRESSURE 







Figure 4 . 107 : Ramp Motion Flow Details. M 0 0 = .3, k = .0127, Re = 2.7 x 10 6 . 
n = 17.00°. 
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DENSITY PRESSURE 






igure 4 . 108 : Ramp Motion Flow Details, M 0 0 = .3, k = .0127, Re = 2.7 x 10*'. 

= 18.00°. 
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DENSITY PRESSURE 




Figure 4.109: Ramp Motion Flow Details. = .3, k = .0127. Re — 2.7 x lOR 
n = 19.00°. 
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DENSITY PRESSURE 




1 igure 4.110: Ramp Motion Flow Details. — .3, k = .0127, Re = 2.7 x 10*'. 
. = '> 0 . 00 °. 
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C. OSCILLATORY MOTION SOLUTION 



1 he unsteady solution for a periodic oscillatory motion, given by a(t) = 4.S6 t- 
J I I at = 0.6. Re c = 4.S x 10'’. with a reduced frequency of k = 0.1(> 

A',h also obtained. Here the reduced frequency is defined as k = ^c/U^. The flow foi 
Mu- motion is initially purely subsonic; bill. <is the angle of attack increases to aboui 
mi / 1 ^5°. supersonic flow conditions are encountered at the leading edge region and a 
'lan-onic shock forms. This shock is present during the upstroke until the maximum 
angle of attack is reached and during the downstroke up to about a(t) ~ 5.0°. The 
ei. m put eel and measured lift and pitching moment response are compared in Figs 
Fill and 4. 112. respectively. 



a(t)=4.86+2.44sin(wt), M=.6, k=.16. Re=4800000 




Angle of Attack, deg 

Figure 4 . 111 : Comparison of Measured and Computed Lift for Oscillatory Test 
Mel ion 

The computed lift and pitching moment coefficients are in close agreement with 
t lie measured values. The computed surface pressure distribution is compared with 
the measurements of reference [6] for two angles during the upstroke and two angle- 
during the downstroke in Figs. 4.113 through 4.116 The computed surface pressure i- 
in better agreement with the measurements at the lower angles of incidence (a = 5.95 



\\> and a = 5.11° down). At higher incidences (Fig. 4.114 and 4.115), the agreement 
'ia! erio rates in the region around the shock. The global view of the computed density 
held shows that the density contours smoothly cross the zonal interface for the case 

ne a shock exists. 



a(t)=4.86+2.44sin(wt), M=0.60, k=0.16, Re=4800000 




Figure 4.112: Comparison of Measured and Computed Moment Coefficient for 

()>< illatorv Test Case 



a=5.95 deg. up, M=0.60, k=0.16, Re=4800000 




Figure 4.113: Comparison of the Measured and Computed Unsteady Surface Pro 
•'lire Coefficient of Oscillatory Motion, o = 5. 95° upstroke. 



a=6.97 deg. up, M=0.60, k=0.16, Re=4800000 




Figure 4.114: Comparison of the Measured and Computed Unsteady Surface Pn 
"lire Coefficient of Oscillatory Motion, a = 6.97° upstroke. 
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a=6.57 deg. down, M=0.60, k=0.16, Re=4800000 




Figure 4.115: Comparison of the Measured and Computed Unsteady Surface Pre- 
-ure Coefficient of Oscillatory Motion, o = fi.oT 0 downstroke. 



a=5.1 1 deg. down, M=0.60, k=0.16, Re=4800000 




Figure 4.116: Comparison of the Measured and Computed Unsteady Surface Pres 
-ure Coefficient of Oscillatory Motion, o = 5.9.")° downstroke. 



V. CONCLUSIONS 



A Mention procedure suitable for steady and unsteady compressible flow sol u - 
tu>n> using zonal overlapped grids was developed. Simple weighted averaging wa- 
nted at the overlapped zonal interfaces. Steady and unsteady, inviscid and visual- 
lluw solutions for subsonic and transonic flows over airfoils were presented to validate 
Hu' zonal grid approach. 

1 he i 11 viscid solutions presented show the overlapped zonal interface's ability tn 
pa-s flow properties without distortion. It was found that for the inviscid test case** 
i he location of the zonal interface is not important. In fact, as the zonal interface 
moved closer to the airfoil, while keeping the number of grid points constant, the 
pressure coefficient prediction actually improved. This is due to more grid point- 
being clustered close to the airfoil. This test case had strong shocks on the upper ami 
M.wer surface of the airfoil, and as the grid points were moved closer to the airfoil, 
the computed solution tended to have bigger oscillations near the shock. 

1 he steady viscous test cases showed that using the Baldwin-Lomax turbulnn <• 
model with the present approach gave accurate results for attached and mildly >cp 
mat i'd flow over stationary airfoils. This case demonstrates one of the advantage- 
wf the present approach, specifically, that solving the inviscid equations on the outer 
m id and the viscous equations on the inner grid gives good results. The flow variable- 
were also passed smoothly through the zonal interface. 

fhe ramp case again showed good agreement with experimental data. With 
this case another advantage of the present approach was displayed. The inner grid 
wih rotated with the airfoil to the new angles of attack. The high order accurate 
m heme enabled the software to convect vortices, and at high angles, these vortices 



convected up to 3 chord lengths. 

The final test case was an oscillatory one. For this test case, supersonic flow was 
encountered on the upstroke at the leading edge region which produced a transonic 
shock. The computed lift and pitching moment coefficients are in close agreement 
with the measured values. At the higher angles of attack the agreement with measured 
data deteriorated in regions near the shock. 

The following are some recommendations based on this study. 

1. It has been shown that the zonal grid approach can be used to study unsteady 
viscous flows. A systematic comparison with other codes should be conducted 
in order to quantify the efficiency of the present approach. 

2. Presently, there exists no easy way to generate zonal grids. In this study the 
grids were generated separately and then refined through several other pro- 
grams. The development of a new or the modification of an existing software 
package such as GRAPE is recommended. The user should be able to specify, 
along with all the usual information, airfoil shape, the location of the overlapped 
region, the number of overlapped cells and the size of the outer grid. 

3. An advantage of the present approach is that it can be extended to multiple 
inner and outer grids. A flow solver needs to developed to take advantage of 
this so that effects of oscillating bodies in relative motion can be studied. 
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